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Abstract 


A fast implicit upwind algorithm for the solution of the time- 
dependent Euler equations is presented for aerodynamic analysis 
involving unstructured dynamic meshes. The spatial discretiza- 
tion of the scheme is based on the upwind approach of Roe re- 
ferred to as flux -difference splitting (FDS). The FDS approach is 
naturally dissipative and captures shock waves and contact dis- 
continuities sharply. The temporal discretization of the scheme 
involves an implicit time- integral ion using a two-sweep Gauss- 
Seidel relaxation procedure. The procedure is computationally 
efficient for either steady or unsteady flow problems. The paper 
gives a detailed description of the implicit upwind solution algo- 
rithm along with results which assess the capability. The results 
are presented for the NACA 0012 airfoil and for the Boeing 747 
aircraft. The 747 geometry includes the fuselage, wing, hori- 
zontal and vertical tails, under-wing pylons, and flow-through 
engine nacelles. Euler solutions for the 747 aircraft on an un- 
structured tetrahedral mesh containing approximately 100,000 
cells were obtained to engineering accuracy in less than one 
hour CPU lime on a Cray-2 computer. 


dynamic analysis involving unstructured dynamic meshes. The 
spatial discretization of the scheme is based on the upwind ap- 
proach of Roe 10 referred to as flux-difference-splitting (FDS). 
The FDS approach is naturally dissipative and captures shock 
waves and contact discontinuities sharply. The temporal dis- 
cretization of the scheme involves an implicit time-integration 
using a two-sweep Gauss-Seidel relaxation procedure. The pro- 
cedure is computationally efficient for either steady or unsteady 
flow problems. The paper gives a detailed description of the im- 
plicit upwind solution algorithm along with results which assess 
the capability. The results are presented for the NACA 0012 
airfoil and for the Boeing 747 aircraft. 

Euler Equations 


In the present study the flow is assumed to be governed 
by the three-dimensional time- dependent Euler equations which 
may be written in integral form as 

J QdV + j ( En r + Fn„ + Gn,)dS = 0 (1) 


Introduction 


In recent years significant progress has been made on de- 
veloping numerical algorithms for the solution of the govern- 
ing fluid flow equations based on unstructured meshes. 1 ITiis 
progress includes improvements in solution accuracy as well as 
computational efficiency. For example, upwind methods have 
been developed for unstructured meshes which are based on 
the local wave propagation characteristics of the flow and con- 
sequently produce highly accurate solutions. 2,3 Most of these 
upwind methods, however, use explicit time-marching schemes 
to integrate the governing equations in time to steady state The 
explicit approach is computationally efficient when applied to 
meshes that are coarse, but the rate of convergence deteriorates 
significantly when finer meshes are used. For cases where finer 
meshes are used, either a multigrid strategy for convergence ac 
celeralion or an implicit temporal discretization which allows 
large lime steps is required to obtain steady-state solutions in a 
computationally efficient manner. Implicit upwind solution al- 
gorithms for unstructured meshes in two dimensions have been 
reported by the author in Ref. 8. These algorithms are simi- 
lar to the point-implicit scheme of Thareja, et al. 9 although the 
methods of Ref. 8 are fully implicit and not point implicit. The 
purpose of the paper is to report the extension of the implicit 
discretization of Ref. 8 to unstructured meshes in three dimen- 
sions. This new flow solver is a fast implicit upwind algorithm 
for the solution of the time-dependent Euler equations, for aero- 
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where the vector of conserved variables Q and the convective 
fluxes F, F, and G are given by 
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The velocities V, V, and IV are defined by 

V = U - xu V ~ v ~ yt' W = w - (3) 


where x h y ( > and c, are the grid speeds in the x,y, and c 
directions, respectively, and the pressure p is given by the 
equation of state for a perfect gas 
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llie above equations have been nondimensionalized by the 
frecstrcam density p^ and the freestream speed of sound 
Also, the second integral in Eq. (1) is a boundary integral result- 
ing from application of the divergence theorem, and n f , n y , and 
are Cartesian components of the unit normal to the boundary 
surface. 

Spatial Discretization 

The spatial discretization is based on Roe’s flux-difference 
splitting which is herein implemented as a cell-centered scheme 
whereby the flow variables are stored at the centroid of each 
tetrahedron and the control volume is simply the tetrahedron 
itself. Consequently, the spatial discretization involves a flux 
balance where the fluxes across the four faces of a given tetra- 
hedron are summed as 

4 \ 
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where A S is the area of the face. The flux vector II is 
approximated by 


distance-weighted average of the flow variables in the tetrahedra 
surrounding node i. The upwind-biased interpolation for q + is 
determined similarly. Also the parameter k in Eq. (7) controls 
a family of difference schemes by appropriately weighting A_ 
and A + . On structured meshes it is easy to show that k = -1 
yields a fully upwind scheme, k = 0 yields Fromm’s scheme, 
and k = 1 yields central differencing. 

On highly stretched meshes, the formula for A+ is modified 
to be 

A + = -T !(«-«>) (9) 

where a and b are the distances from the midpoint of the face to 
the centroids of tetrahedra j and respectively. This formula 
weights the flow variables in the interpolation formula (Eq. 
(7)) differently to account for the stretching of the mesh. For 
example, by substituting Eq. (9) into Eq. (7) and letting k = 1 
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Furthermore, in calculations involving upwind-biased schemes, 
oscillations in the solution near shock waves are expected to 
occur. To eliminate these oscillations flux limiting is usually 
required. The flux limiter modifies the upwind-biased interpo- 
lations for and </ + such that, for example 


q~ =(ij F - ^)A. F (I F *«$)A+] (11) 


where .s is the flux limiter. In the present study, a continuously 
differentiable flux limiter was employed which is defined by 
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^ where ( is a very small number used to prevent division by zero 
in smooth regions of the flow. 


where Q~ and Q ¥ are the state variables to the left and right 
of the cell face and A is the flux Jacobian matrix given by 
<)II/()Q. Also the tilde and the absolute value sign indicate that 
the flux Jacobian is evaluated using the so-called Roe-averaged 
flow variables and the absolute value of the characteristic speeds. 


The left and right states Q and Q * , are determined by 
upwind-biased interpolations of the primitive variables </. In 
three dimensions, for a given tetrahedron j, for example, the 
upwind-biased interpolation for across the common face 
between tetrahedra j and k is defined by 

</ ' ~ <tj F -[(I - k)A_ F ( l F k)A j.) (7) 

where 

A+ - r/jL - <lj 

A - - Hj - f/« (8b) 

In Eqs. (7) and (8), i\ } and < ik arc the vectors of primitive 
variables at centroids j and respectively, and #/„ the vector of 
primitive variables at node i (the node of tetrahedron j opposite 
to the face being considered), is determined by an inverse- 


Temporal Discretization 

The temporal discretization is an implicit time-marching 
scheme involving a Gauss-Seidel relaxation procedure. The 
scheme is derived in general by first linearizing the flux vector 
II according to 

//” +1 = //" + ^A Q (13) 


where ()II/DQ is the flux jacobian A, as discussed before, and 
AQ — Q tt+I — Q n . Linearizing both flux terms on the right- 
hand-side of Eq. (6) using Eq. (13), and ignoring the tilde on 
the flux jacobian, results in 
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where / is the identity matrix, “vol” is the volume of the tetra- 
hedron j, and A Q m is the change in flow variables in each 
of the four tetrahedra adjacent to tetrahedron j. Also in Eq. 
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(14), A + and A~ arc forward and backward llux jacobians, re- 
spectively. For flux-difference splitting, the exact Jacobian A 
(derivative of the right-hand-side of Eq. (6) with respect to Q) 
is too expensive to compute and thus an approximate Jacobian 
is normally used. This is accomplished by constructing the ja- 
cobians making use of the fact that the forward and backward 
jacobians should have non-negative and non-positive eigenval- 
ues (characteristic speeds), respectively. This is accomplished 
by expressing alternatively the jacobians using similarity trans- 
formations such that 

A + - AM 4 /r’ . A~ - WA‘/r’ (15) 

where A 1 and A" are diagonal matrices whose diagonal ele- 
ments are the eigenvalues A + and \~ defined by 

A + = 1(A + | A | ) A _ =’(A-|A|) (16) 

and H is ihe matrix whose columns are the corresponding 
eigenvectors. 

Direct solution of the system of simultaneous equations 
which results from application of Eq. (14) for all tetrahedra 
in the mesh, requires the inversion of a large matrix with 
large bandwidth which is computationally expensive. Instead, a 
Gauss-Scidel relaxation approach is used to solve the equations 
whereby the summation involving A Q,„ is moved to the right 
hand side of Eq. (14). The terms in this summation are then 
evaluated for a given time step using the most recently computed 
values for AQ,„. The solution procedure then involves only the 
inversion of a 5 x 5 matrix (represented by the terms in square 
brackets on the left hand side of Eq. (T4)) for each tetrahedron 
in the mesh. Also, although the procedure is implemented for 
application on (randomly-ordered) unstructured meshes, it is not 
a point Gauss-Seidel procedure. The method is in fact more like 
line Gauss-Seidel since the list of tetrahedra that make up the 
unstructured mesh is re-ordered from upstream to downstream, 
and the solution is obtained by sweeping two times through the 
mesh as dictated by stability considerations The first sweep is 
performed in the direction from upstream to downstream and 
the second sweep is from downstream to upstream. For purely 
supersonic flows the second sweep is unnecessary. ^ 


Boundary Conditions 


To impose the (low tangcncy boundary conditions alo ng the 
surface of the vehicle, the flow variables are set within dummy 
cells that arc effectively inside the geometry being considered. 
The velocity components within a dummy cell, (u,i',u>)j. are 
determined from the values in the cell j adjacent to th e su rface, 

t„ . This is accomplished by first rotating the components 

into a coordinate system that has a coordinate direction normal 
,« the boundary face. The sign of the velocity component in 
this direction is changed (hence imposing no flow through the 
face) and the three velocity components are then rotated back 
into the original .r.j coordinate system. Af.cr considerable 
algebra this yields 
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where n x , n„, and n, are the x, y, and « components of the unit 
vector that is normal to the boundary face. Also, pressure and 
density within the dummy cell are set equal to the values in the 
cell adjacent to the surface. 

After application of the upwind-biased interpolation formula 
to determine and ^ at each face, the velocity components 
are corrected to give a "strong” implementation of the surface 
boundary condition according to 

n«,r«rfrj = « “ n r(» n x + ” U !I + «’"*) 

Vcorrtct'd = “ n y (un, + vn y + urn,) (18) 

WcoTTtiUd = - n z {un z + vn y -f wn,) 

In the farfield a characteristic analysis based on Riemann 
invariants is used to determine the values of the flow variables on 
the outer boundary of the gnd. This analysis correctly accounts 
for wave propagation in the farfield which is important for 
rapid convergence to steady-state and serves as a “nonrcfleciing 
boundary condition for unsteady applications. 

Results and Discussion 

To assess the accuracy and efficiency of the implicit up- 
wind solution algorithm, calculations were first performed in 
two dimensions for the NACA 0012 airfoil. These results were 

obtained using the unstructured mesh shown in Fig. 1. The grid 

has 3300 nodes and 6466 triangles, and extends 20 chordlengths 
from the airfoil with a circular outer boundary. Also there are 
110 points that lie on the airfoil surface. Steady-state calcula- 
tions were performed for the airfoil at a freestream Mach number 
of Moo = 0.8 and an angle of attack of a = 1.25°. The results 
were obtained using both the implicit relaxation time-marching 
scheme and an explicit three-stage Runge-Kutta time-marching 
scheme. The explicit results were obtained using a CFL number 
of 4.0, with residual smoothing and local time-stepping to ac- 
celerate convergence to steady state. The implicit results were 
obtained using a CFL number of infinity. Such a large value 
was used for the implicit results since the relaxation scheme 
has maximum damping and hence fastest convergence for very 
large time steps. This is in contrast with implicit approximate 
factorization schemes which have maximum damping for CFL 
numbers on the order of 10. Also, the flux jacobians of the 
implicit scheme were updated only every twenty iterations. 

A comparison of the convergence histories between explicit 
and implicit time-marching is shown in Fig. 2(a). The error 
in the solution was taken to be the L 2 norm of the density 
residual. As shown in Fig. 2(a), the explicit solution is slower to 
converge than the implicit solution. The explicit solution takes 
approximately 739 CPU secs. (2,682 iterations) on a Cray-2 
computer to converge to engineering accuracy, which is taken 
to be a three order of magnitude reduction in solution error. In 
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Fig. I Partial view of unstructured mesh of triangles about 
the NACA 0012 airfoil. 


constrast, the implicit solution is converged to three orders of 
magnitude in only approximately 362 secs. (1,251 iterations). 
The resulting steady pressure distribution is shown in Fig. 2(b). 
For this case there is a relatively strong shock wave on the 
upper surface of the airfoil near 62% chord and a relatively 
weak shock wave on the lower surface near 30% chord. The 
pressure distribution indicates that there is only one grid point 
within the shock structure, on either the upper or lower surface 
of the airfoil, due to the sharp shock capturing ability of the Roe 
solver. 

To assess the efficiency of the implicit upwind solution 
algorithm in three dimensions, calculations were performed for 
the Boeing 747 aircraft. These results were obtained using the 
unstructured mesh shown in Fig. 3. The 747 geometry includes 
the fuselage, the wing, horizontal and vertical tails, under-wing 
pylons, and flow-through engine nacelles. The unstructured 
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(a) convergence histories. 
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(b) pressure distribution. 

Fig. 2 Comparison of steady-state results for the NACA 
0012 airfoil at Moo = 0.8 and o = 1.25°, 
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whereas the implicit solution required less than half of that, or 
3,578 secs, (522 steps). The resulting steady pressure coefficient 
contours on the surface of the 747 aircraft are shown in Fig. 
5. The contours indicate that there is a significant amount of 
flow compression on the nose of the aircraft, along the inboard 
leading edge of the wing, and inside the cowl of the engine 
nacelle. There is flow expansion on the forward fuselage, on the 
horizontal and vertical tail surfaces, and on the upper surface of 
the wing terminated by a shock wave. Additional details of the 
mesh and pressure contours on the outboard pylon and engine 
nacelle are shown in Fig. 6. These contours show further flow 
expansion on the outside of the cowl and within the inner core 
of the engine. 


Fig. 4 Comparison of convergence histories for the Boeing 
747 aircraft at U* = 0.84 and a = 2.73°. 


Concluding Remarks 


mesh for the 747 contains 101,475 tetrahedra and 19,055 nodes 
for the half-span airplane. Also there are 4,159 nodes and 
8,330 triangles on the boundaries of the mesh which include 
the airplane, the symmetry plane, and the farfield. Steady-state 
calculations were performed for the aircraft at 3/^ = 0.84 and 
o = 2.73°. The results were obtained using both the implicit 
and explicit lime-marching schemes. Similar to the NACA 0012 
cases, the explicit results were obtained for the 747 using a 
CFL number of 4.0, with residual smoothing and local time- 
stepping to accelerate convergence to steady state. The implicit 
results were obtained using a CFL number of infinity and the 
flux jacobians were updated only every twenty iterations. 

A comparison of the convergence histories between explicit 
and implicit time-marching is shown in big. 4. The explicit 
solution required 8,251 CPU secs. (1,124 iterations) on a Cray- 
2 computer to converge the solution three orders of magnitude, 


A fast implicit upwind algorithm for the solution of the 
lime-dependent Euler equations was presented for aerodynamic 
analysis involving unstructured dynamic meshes. The spatial 
discretization of the scheme is based on the upwind approach 
of Roe referred to as flux-difference splitting (FDS). The FDS 
approach is naturally dissipative and captures shock waves and 
contact discontinuities sharply. The temporal discretization of 
the scheme involves an implicit time- integration using a two- 
sweep Gauss-Seidel relaxation procedure. The procedure is 
computationally efficient for either steady or unsteady flow prob- 
lems. Results were presented for the NACA 0012 airfoil and 
for the Boeing 747 aircraft. The 747 geometry included the 
fuselage, wing* horizontal and vertical tails, under- wing pylons, 
and flow-through engine nacelles. Euler solutions for the 747 
aircraft on an unstructured tetrahedral mesh containing approx- 
imately 100,000 cells were obtained to engineering accuracy in 
less than one hour CPU time on a Cray-2 computer. 



Fig. 5 Steady pressure coefficient contours on the Boeing 747 aircraft at 
M^ t - 0.84 and a = 2.73°. 
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